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Simulating local adaptation to climate of forest trees with a 
Physio-Demo-Genetics model 



One challenge of evolutionary ecology is to predict the rate and mechanisms of 
population adaptation to environmental variations. The variations in most life 
history traits are shaped both by individual genotypic and by environmental vari- 
ation. Forest trees exhibit high levels of genetic diversity, large population sizes, 
and gene flow, and they also show a high level of plasticity for life history traits. 
We developed a new Physio-Demo-Genetics model (denoted PDG) coupling (i) 
a physiological module simulating individual tree responses to the environment; 
(ii) a demographic module simulating tree survival, reproduction, and pollen 
and seed dispersal; and (iii) a quantitative genetics module controlling the herita- 
bility of key life history traits. We used this model to investigate the plastic and 
genetic components of the variations in the timing of budburst (TBB) along an 
elevational gradient of Fagus sylvatica (the European beech). We used a repeated 
5 years climatic sequence to show that five generations of natural selection were 
sufficient to develop nonmonotonic genetic differentiation in the TBB along the 
local climatic gradient but also that plastic variation among different elevations 
and years was higher than genetic variation. PDG complements theoretical mod- 
els and provides testable predictions to understand the adaptive potential of tree 
populations. 
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Introduction 

The ongoing and predicted rapid changes in temperature, 
precipitation and CO2 atmospheric concentration, and the 
resulting increase in the frequency and intensity of extreme 
events such as droughts will have a wide range of long-term 
implications for natural population dynamics and ecosys- 
tem sustainability. Within a population, these changes 
impose strong selective pressures, which affect demo- 
graphic rates and can cause genetic evolution across gener- 
ations (Hansen et al. 2012). Moreover, climate change 
(CC) also affects the physiology and development of indi- 
vidual organisms up to the limits of their phenotypic plas- 
ticity, that is, the ability of individual genotypes to produce 
alternative phenotypes in different environments (Chevin 
et al. 2013). The interplay between genetic evolution and 
phenotypic plasticity ultimately determines a population's 
ability to adjust (without migrating) to novel environmen- 
tal conditions imposed by CC. Investigating these mecha- 
nisms is essential for predicting eco- evolutionary dynamics 



and ecosystem services and for guiding conservation 
efforts. 

This issue is crucial for trees because of their pivotal role 
in the functioning and biodiversity of forest ecosystems. 
Multisite experiments (using forester provenance tests) 
showed that current tree populations can adjust to varying 
environmental conditions through phenotypic plasticity 
over a non-negligible latitudinal range (Rehfeldt et al. 
2002). Moreover, based on the patterns of local adaptation 
displayed by most tree species over the course of postglacial 
recolonization, forest tree populations are usually assumed 
to have a high evolutionary potential (Savolainen et al. 
2007; Alberto et al. 2013); however, tree population abili- 
ties of genetic evolution over a short timescale (i.e., micro- 
evolution) remain largely unresolved. In addition, plasticity 
and genetic adaptation can interact together and with gene 
flow, as illustrated mostly by theoretical models or studies 
of model species. One well-known interaction is the inter- 
play between gene flow and adaptation when the environ- 
ment changes both in space and in time (Polechova et al. 
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2009). Trees are capable of long-distance pollen-mediated 
gene flow, which could promote adaptive evolution to 
novel environments (Kremer et al. 2012). Another perva- 
sive interaction involves plasticity and genetic adaptation; 
plasticity can be adaptive if plastic trait variation increases 
individual fitness (Nicotra et al. 2010), or it can be mal- 
adaptive if plasticity decreases fitness (Ghalambor et al. 
2007). Moreover, when adaptive plasticity cannot evolve, it 
can slow down the genetic response to directional selection, 
but it also allows phenotypes to track environmental 
change more closely (Chevin et al. 2013). 

Methodological developments currently limit our under- 
standing of the interplay among plasticity, genetic adapta- 
tion, and gene flow and their impacts on tree population 
dynamics. In most evolutionary models thus far developed 
for tree life history traits, individual fitness is either directly 
controlled by the genotype (Le Corre and Kremer 2003) or 
derived from genetically controlled life history traits (Ku- 
parinen et al. 2010). Climate effects on water and carbon 
exchanges are a complex process that has been studied by 
ecophysiologists and have rarely been explicitly taken into 
account as a selective pressure in evolutionary models (but 
see (Kramer et al. 2008). Similarly, the interindividual vari- 
ation and adaptive potential of traits related to climate 
response have rarely been incorporated into biophysical 
and ecophysiological models. 

The different time scales considered by ecophysiological 
and evolutionary models (typically from the hour to the 
year or tens of years for the former versus many genera- 
tions at equilibrium for the latter) are generally considered 
to be challenges in the development of coupled physio- 
demo-genetic models. However, neither of these time 
scales may be relevant for forming accurate predictions of 
realistic tree population responses to CC. Indeed, current 
forest tree populations can rarely be considered to be at 
equilibrium, and demographic processes play a major role 
in the dynamics of adaption over a few generations (Savo- 
lainen et al. 2007). Moreover, CC is likely to involve non- 
continuous and nonpredictable change in response to 
abiotic conditions, which limit the relevance of long-term 
predictions at equilibrium (Kremer and Le Corre 2011). 
In contrast, the predictions of biophysical and ecophysio- 
logical models cannot be generalized over more than one 
generation if the microevolution of functional traits within 
a population is not negligible. Therefore, there is a need 
for physio-demo-genetic models to address the timescale 
of a few generations (<10). 

In this study, we propose a new Physio-Demo-Genetics 
model (PDG) coupling the following: (i) a functional mod- 
ule derived from CASTANEA (Dufrene et al. 2005) to sim- 
ulate carbohydrate and water fluxes at the tree level using 
daily climate observations; (ii) a population dynamics 
module to convert carbohydrate reserves into demographic 
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rates for adult trees (growth, mortality, and seed produc- 
tion) and to simulate ecological processes across the life 
cycle (including seed and pollen dispersal, germination 
rate, and density-dependent mortality of seedlings); and 
(iii) a quantitative genetics module relating genotype of the 
quantitative trait loci (QTL) to the phenotype of one or 
more functional traits (Labonne and Hendry 2010). This 
individual-based, spatially explicit model simulates the evo- 
lution of functional traits in tree populations, where phe- 
notypic differences between individuals are determined by 
their genotype at QTLs that control functional traits and by 
their physiological response to local climate conditions. 
PDG is available from the CAPSIS modeling platform (Du- 
four-Kowalski et al. 2012). 

We used PDG to simulate local adaptation in a continu- 
ous tree population that expands along an elevational gra- 
dient based on experimental data collected from Fagus 
sylvatica populations on Mont Ventoux in southeastern 
France. Under divergent selection, local adaptation is 
expected to result in phenotypic differentiation for traits 
contributing to fitness across the gradient. Among the vari- 
ous traits contributing to fitness, we focused on the timing 
of budburst (TBB), a phenological trait that determines the 
length of the growing season in F. sylvatica (Davi et al. 
2011). An earlier budburst extends the time during which 
photosynthesis occurs (Richardson et al. 2006), but it also 
increases the risk of late frost damage (Dittmar and EUing 
2006). In situ observations of F. sylvatica on Mont Ventoux 
revealed a classical phenotypic cline in TBB resulting from 
plastic variation; budburst occurs earlier at lower than 
higher elevations because TBB is triggered by the heat sum 
(Davi et al. 2011). The opposite is observed for genetic 
clines as assessed in common garden experiments, where 
TBB is observed under the same environmental conditions 
for all populations; under such conditions, _F. sylvatica 
populations originating from higher elevations show earlier 
budburst than those originating from lower elevations (von 
Wuehlisch et al. 1995; Vitasse et al. 2009a; Gomory and 
Paule 2011). This situation in which the phenotypic and 
genetic clines vary in opposite directions is referred to as a 
counter-gradient variation. In contrast, genetic and pheno- 
typic clines have been shown to exhibit co-gradient varia- 
tion for TBB in some species (e.g., Quercus sp.), while clear 
linear genetic clines are not observed for other species (e.g., 
Fraxinus) (Vitasse et al. 2009a). 

In this article, we illustrate the potential of PDG to eluci- 
date the processes through which adaptation proceeds for 
_F. sylvatica on Mont Ventoux. We address the following 
issues: (i) How do adaptive genetic variation and pheno- 
typic plasticity contribute to TBB variation along an eleva- 
tional gradient? (ii) How fast can genetic differentiation in 
TBB develop? and (iii) Is there a monotonic trend in the 
genetic variation across the gradient? 
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Materials and methods 

Overview of PDG model 

The physiological process-based module 
This module corresponds to tlie CASTANEA library hosted 
on the CAPSIS platform. Initially developed at the stand 
scale, CASTANEA (Dufrene et al. 2005) simulates canopy 
photosynthesis (i.e., gross primary production, GPP) and 
transpiration, maintenance and growth respiration, sea- 
sonal development and assimilate partitioning to leaves, 
carbohydrate storage (hereafter reserves), stems, branches, 
and coarse and fine roots. The meteorological driving vari- 
ables are global radiation, rainfall, wind speed, air humid- 
ity, and temperature. A complete description of the model 
is given in Dufrene et al. (2005), and the submodel of car- 
bon allocation is described by Davi et al. (2009). 

In its initial version, CASTANEA simulated CO2 and 
H2O fluxes, considering one average tree as representative 
of the whole stand. To account for interindividual varia- 
tion, we considered each tree as a single unit with its own 
parameters for the CASTANEA simulation. Note that all 
CASTANEA units were treated independently, meaning 
that we do not account for competition among trees for 
light or soil water acquisition. In contrast to the stand- 
level version, we computed several variables at the indi- 
vidual tree level, including biomass (5tree)' leaf area index 
(LAI) by tree, and crown projection {A^^avin)^ to deter- 
mine the carbon and water budgets of each tree (see 
Appendix SI). 

Only one of all the CASTANEA parameters was allowed 
to vary among trees, namely, f critBB' the critical value of the 
state of forcing, which is most commonly referred to as the 
temperature sum required for budburst. The TBB was sim- 
ulated following the equs 1-3: 

_ / if r > and N > Nstartl 1 

WfrcBB - I 0 if r < 72 or N < Nstartl J 

(1) 

where RtcBB is the rate of forcing for bud break, T the mean 
daily temperature, T2 the base temperature, N the day of 
year, and Nsxarti the date of onset of rest. 

N 

SfrcBB = J^frcBB if SfrcBB < fcritBB (2) 

NsTAfiTl 

TBB = NifSfaBB>FcritBB (3) 

with Sfi-cBB the state of forcing, -FcritsB the critical value of 
state of forcing for the transition from quiescence to the 
active period and TBB the day when bud break occurred. 

We chose to focus on the fcritBB variation for two rea- 
sons. First, -FcritBB is one of 17 key parameters controling 



carbon exchange according to the sensitivity analyses per- 
formed in Dufrene et al. (2005). Second, the genetic clines 
observed for TBB in a common garden of populations 
from different elevations are likely to result from among- 
population variation for -FcritBB- For the sake of simplicity, 
we assumed that the effect of chilling on dormancy was 
constant across year and elevations. 

To assess the impact of drought on photosynthesis, an 
annual water stress index (WSI) was estimated at the tree 
level as the sum of the daily reduction of photosynthesis 
caused by soil drought (i.e., when the relative soil water 
content drops below 40% of the soil water reserve). 

In this version of CASTANEA, we also added the effect 
of late frost on LAI. Trees are most sensitive to frost during 
budburst; frost can kill the new shoots, reduce growth, and 
cause misshapen branching (Dittmar and EUing 2006). 
When unfolding leaves are affected by frost, decreased leaf 
area is expected. In PDG, every day during which the mini- 
mal daily temperature (T^in) fell below a threshold value 
(rminEffect) foUowing the initiation of budburst was consid- 
ered to affect the LAI as follows: 

LAItree(day-F 1) = LAItee(day) 

(1 + (T'min Effect — rmm)CoefFrostE£fect) 

Parameterization. In previous studies, species-specific CAS- 
TANEA parameters for F. sylvatica were determined and 
CASTANEA was validated at an experimental site in 
Northern France (Davi et al. 2005). Additionally, some 
site-specific parameters were measured in Ventotix 
(Table 1). First, the budburst model was calibrated using 
two types of data. The onset date of rest (-RfrcSB) was esti- 
mated by an experiment on dormancy release in the spring 
of 2012 (unpublished data). The average critical value for 
the state of forcing (.FcritBB) was estimated using budburst 
survey from 2007 to 2011 at two elevations (1117 and 
1340 m) for 20 trees per elevation (Fig. 1). Characteristics 
of sun leaves (nitrogen content and LMA = leaf mass per 
area) were obtained for 149 trees in an intensive-studied 
plot (Bontemps 2012). Canopy clumping (CI) was esti- 
mated using five hemispherical photographs taken in 
the same plot in the summer of 2008, foUowing the 
methodology described in Davi et al. (2008). The 
photosynthesis parameter (maximum carboxylation 
rate = 51.6 |imol photon m^^ s^^) was estimated from 
previous measurements at the same site in the summer of 
2006 (M. Ducrey and R. Hue, personal communications). 

The demographic module 

Adult growth, mortality, and seed production. The reserves 
produced by photosynthesis at a daily time step were allo- 
cated to growth and used to predict tree mortality. Two 



©201 4 The Authors. Evolutionary Applications pubWshad, by John Wiley & Sons Ltd 7 (2014) 453-467 



455 



Physio-Demo-Genetic model for trees 
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Parameter 


Acronym 


Value 


Unit 


Sources 


Pliysical/pliysiological CASTANEA module 










Canopy clumping 




0.56 


- 


Davi et al. (2008) 


Nitrogen content 




2.2 


% 


Bontemps (2012) 


Leaf mass per area 




93 




Bontemps (2012) 


Relationship between maximal rate of carboxylation 




26.04 


Ijmol CO2 g N"' s"^ 


Ducrey and Hue (personal 


and nitrogen 








communications) 


Date of rest onset for budburst 




78 


Days 


Davi (unpublished data) 


Average critical value of forcing state 




190 


°C 


This study 


Base temperature for forcing budburst 




0 


°c 


Fixed 


Ratio between fine root and leaf biomass 




1 


- 


Fixed 


Soil extractable water 


SEW 


60 


mm 


Nourtieretal. (2012) 


Threshold value for frost effect on leaf area index 


^minEffect 


0 


"C 


Fixed 


Demographic module 










Critical threshold of carbon reserves at the end of the 


sSres 


100 


gC msoif^ 


Fixed 


year for reproduction 










Critical threshold of carbon reserves at the end of the year 


CumCR 


45 


gC msoif^ 


Fixed 


Maximal difference between carbon needs and 


bbCR 


160 


gC msoif^ 


Fixed 


carbon reserves before budburst 










Rate of seed production 


/?SP 


0.05 


- 


Fixed 


Cost to produce one seed 


C 


0.45 


gc 


Han etal. (2011) 


Rate of empty seeds 


fES 


0.33 


- 


Oddou-Muratorio (personal 
observations) 


Rate of seed germination 




0.485 


- 


Oddou-Muratorio (personal 
observations) 


Rate of seed survival 


fss 


0.15 


- 




Average distance of seed dispersal 


As 


18.13 


m 


Bontemps et al. (2013) 


Shape of the seed dispersal kernel 


bs 


0.31 




Bontemps et al. (2013) 


/nvtri clutr uiiLdiiLtr ui [JUiitrii uij[Jfcii !>cii 


Op 


37 9 




Rnntpmns pt ^\ ('?f)^^\ 


Shape of the pollen dispersal kernel 


bp 


0.97 




Bontemps et al. (2013) 


Parameter relating tree diameter and male fertility 


I'm 


0.82 




Bontemps et al. (2013) 


Rate of selfing 


s 


0.025 




Bontemps et al. (2013) 


Mean height of newly recruited tree 




9 


m 


P. Dreyfus (personal communications 


Standard deviation of newly recruited tree height 


SDh 


0.34 


m 


P. Dreyfus (personal communications 


Mean diameter at breast height (DBH) of newly 


I^DBH 


13.8 


cm 


P. Dreyfus (personal communications 


recruited tree 










Standard deviation of newly recruited tree DBH 


SDdbh 


0.9 


cm 


P. Dreyfus (personal communications 



levels of carbon reserves were considered: the carbon 
reserve at the end of the year (CumCR; Table 1) and the 
difference between the carbon reserves before budburst and 
the amount of carbon required for the complete develop- 
ment of new leaves (bbCR; Table 1). Below a critical level 
of one of these two indicators, a tree would die. Critical lev- 
els of CumCR and bbCR were estimated based on mortality 
rates assessed on Mount Ventoux (H. Davi unpublished 
results). 

Biomass allocated to wood between the date of budburst 
and leaf senescence was converted into a diameter at breast 
height (DBH) increment. Finally, at the end of the year, if 
the biomass of accumulated reserves (Bres) exceeded the 
critical rate for seed production (si5i.es)) the reserve was 
converted into primary seed production (Ns) for each tree 
as follows: 



3 

T3 
S 

E 

CO 







y = 0.99x + 2.74 




R^ = 0.84 










0 1117 m 




* 1340 m 



80 



90 100 110 120 130 140 

Measured budburst 



Figure 1 Observed budburst in Mont Ventoux versus simulated bud- 
burst using an average value of 190°C for the temperature sum 
required for budburst (fcritBs)' 
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Ns = Cpj^sp( ^"^ (5) 

where Cp is the crown projection of the tree (Appendix 
SI), -Rsp is the rate of seed production (Table 1), and c is 
the cost to produce one seed. Parameter c was estimated 
using the dry mass and carbon content of seeds and cupules 
(Han et al. 2011) assuming an associated respiratory cost 
of 50%. 

The effective seed production of a tree, that is, its female 
fecundity, was computed as follows: 

_F = Ns(l - tes) rssrsG (6) 

where r^s is the rate of empty seeds, r^s is the rate of seed 
survival, and rsc is the rate of seed germination. Parameters 
tes and rsc were calibrated on the basis of a germination 
experiment in the years 2009-2010, during which 60 seed 
lots were collected at three elevational levels from Mont 
Ventoux (20 mother trees per altitude level), and from 100 
to 300 seeds per mother tree were scanned (to measure the 
Tes) and sown after stratification (to measure the rsc). 

Pollen dispersal and mating. Pollen dispersal was modeled 
by an exponential dispersal kernel describing the probabil- 
ity that a pollen grain emitted at position (0,0) would polli- 
nate a seed tree at distance r as follows: 

pp(ap,r) = ^ exp( -— ) (7) 
271 ■ flf \ apy 

where the scale parameter flp is homogeneous to the mean 
distance travelled by pollen grain ((5p) with the relationship 
(5p = 2ap (Table 1). The simulation domain was defined 
with reflecting borders to avoid the loss of border tree 
progeny. 

The contribution Ttjj- of poUen tree k to the outcrossing 
poUen cloud to the fertilization of seed tree j (jj^k) 
depended on the distance between trees j and k through the 
poUen dispersal kernel pp and on its pollen fecundity as 
related to the diameter DBHj^ as follows: 

7r^i = pp(ap,bp) xe-D""* (8) 

where )'„ is related to the diameter effects on male fertility. 
We assumed no pollen limitation. 

In addition to outcrossing, selfing was considered to 
occur at a fixed rate s. Tree j was self-pollinated with proba- 
bility 5 and pollinated by other individuals with probability 
(1— s) Kjp. (for 1 < k<N and k=/^j), where N is the total adult 
population size. Parameters 5p, s, and were estimated 
for the _F. sylvatica trees on Mont Ventoux (Oddou-Mura- 
torio et al. 2010). 

Seed dispersal and recruitment (density-dependent mortal- 
ity). The simulation domain was divided in squared cells 
to model seed dispersal. The intensity of seed rain from 



a given seed tree j on the center of cell i was expressed 
by 

^1) = ps{a„rij)Fj (9) 

where ps is the seed dispersal kernel, r^ is the distance from 
tree to the center of cell and Fj is the female fecundity of 
seed tree j. The seed dispersal kernel was modeled as pollen 
dispersal using the exponential kernel described by eqn 4. 
From Ty, we computed the number of new trees from a 
given seed tree j on the whole cell i of area S, as detailed in 
Appendix S2. 

Within each cell i, ^jN^ individuals (j = 1 to the total 
number of seed trees) were created at the age of 40 years, 
thus assuming that the phenology selection did not proceed 
differently before and after this life stage. The height and 
diameter of newly created trees were drawn in a Gaussian 
distribution of parameter {i^h; SDh} and {|.Idbh; SDobh}. 
respectively. The spatial position of each new tree was allo- 
cated randomly within the cell unless its mother tree was in 
the same cell i; in this case, spatial positions were drawn in 
a Gaussian dispersal kernel around the position of the 
mother tree. Mortality during recruitment was modeled as 
a spatial, random (i.e., independent from TBB), density- 
dependent process considering that no tree pairs could 
have more than 30% overlapping crown. 

A quantitative genetic module for the TBB 
The variation in TBB between individuals depended on 
both (i) the individual genetic variation in fcritBs and (ii) 
environmental variation among individual locations and 
years. In a given environment, the higher is the -FcritBB; the 
later is the TBB. In most of the simulated scenarios, the 
environmental component of TBB variation was fuUy 
determined by the variation in daily temperatures during 
the spring, which varied across years for a given individual, 
and across elevations for different individuals. Note that in 
our simulations, elevation could vary by 200 m between 
individuals within the same population (see paragraph 
'Simulation result analyses' below), and thus, both the vari- 
ations in -FcritBB and elevation contribute to the variation in 
TBB within population. 

The value of -FcritBB was determined by ten independent 
biallelic loci with purely additive effects. The occurrences 
of mutations and new allele immigration from other pop- 
ulations than those simulated were ignored. The contribu- 
tion of a genotype at a given locus 1 to FcritBB was given as 
mi+ai, Mil, and mi— ai for the homozygote AlAl, the het- 
erozygote A1A2, and the homozygote A2A2, respectively. 
All the miS were identical and equal to M-f^^imj/lO (Table 
SI). We followed the method proposed by Bost et al. 
(2001) to generate the distribution of allelic effects. They 
showed that for Ni loci having an equal m; = |j,/Ni contri- 
bution to the trait, an L-distribution of QTL effects could 
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be simulated with allele effects randomly drawn from a 
Gaussian distribution of mean |j,/2Ni and with a standard 
deviation a small enough to ensure that di belongs to [0; 
l-i 1^'] (here, a = |j,/8Ni). Allelic effects cxi were constant 
over time (Table 1). The genetic variance for -FcritBB within 
the population/subpopulation was computed as follows: 

Va = — pi)flp where pi is the frequency of allele 

Al at locus 1. 

In most scenarios with microevolution, we considered 
J^critBB to be fully heritable; the narrow-sense heritability 
/j^5 of -FcritBB was Set to 1, meaning that the phenotypic var- 
iance of -FcritBB (^p) equaled to the additive genetic variance 
Va (as h\^^ = Va/Vp), and that the environmental variance 
Ve was 0. Each individual Fcri,BB value was the sum of 
genotypic contributions (see above) across the 10 loci. In 
the control scenario, FcritBB was variable but not heritable 
(i.e., Va = 0); to match the phenotypic variation obtained 
in microevolution scenarios, Vp was set to Vg = 22. Indi- 
vidual FcritBB values were randomly drawn from a Gaussian 
distribution of mean and of variance Vp. Finally, we 

also considered the case where FcritBB was itself a quantita- 
tive trait with = 0.6 (Kramer et al. 2008); each individ- 
ual FcritBB value was the sum of an additive genetic 



component and of a stochastic environmental component, 
drawn in a Gaussian distribution of mean 0 and of variance 
Ve, so that Ve = (1 - /jL) Vp and Va = 22. 



Simulation design and testing hypothesis 

We applied PDG to simulate the evolutionary dynamic of 
F. sylvatica along an elevational gradient from 700 to 
1700 m on a 20 ha grid (200 x 1000 m) divided into 500 
cells (20 X 20 m) (Fig. 2). This case study mimic an eleva- 
tional gradient located on Mont Ventoux (southeastern 
France, 44°10'28"N; 5°16'16"E), where F. sylvatica recently 
expanded under the black pines (Pinus nigra) that were 
planted at the end of the nineteenth century. The species 
currently extend from 750 to 1700 m in elevation on the 
northern aspect. Environmental, climatic and ecological 
data were available from previous studies (Cailleret and 
Davi 2011; Davi et al. 20II). 

Initial conditions and neutral pre-evolution 
The initial population included 500 trees (all 40 years old). 
Their initial height and diameter were drawn by following a 
Gaussian distribution with parameters {[tn; SDh} and 
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Figure 2 Spatial population dynamics over five generations (from GO to GB) under 'adaptive evolution' scenario (B) along the elevational gradient. 
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{Hdbh; SDdbh} (Table 1). Initial allele frequencies were 
drawn from a uniform distribution [0,1], and they formed 
genotypes according to the Hardy-Weinberg equilibrium. 
The initial population was split into five equal subpopula- 
tions at the following five elevational ranges: 790-810, 990- 
1010, 1190-1210, 1390-1410, and 1590-1610 m (Fig. 2). 
To introduce initial genetic differentiation among subpop- 
ulations and spatial genetic structure within subpopulations 
similar to the levels observed on the field, we first simulated 
five generations without selection in which the allelic fre- 
quencies within and among subpopulations evolved only as 
a result of genetic drift and gene flow (Appendix S3). 

Adaptive evolution over six generations 
After initialization, we simulated six nonoverlapping gener- 
ations (from GO to G5) of evolution with selection in addi- 
tion to genetic drift and gene flow. Selection occurred 
within each generation both through differential mortality 
and differential reproduction of individual trees. Each gen- 
eration lasted 70 years and included (i) a seedling stage 
(from years/ages 0 to 39) and (ii) an adult stage (from years/ 
ages 40 to 70). Adult trees grew for 20 years without repro- 
ducing (i.e., until age 60) and then grew and reproduced 
over the course of 10 years. After 70 years, all the surviving 
trees were removed. Seeds produced during generation Gx 
were sent into dormancy until the beginning of generation 
Gx+i, when survival depended on the total seedling density. 

Climatic data 

The basis climate variables (Xbasis) were assessed from 2002 
to 2006 using daily meteorological data measured at a per- 
manent weather station located in Ventoux (Porte et al. 
2004). The elevation effects on temperature, relative 
humidity and precipitation were estimated using Hnear 
models and data acquired from April 2007 to October 2009 
using five HOBO Pro V2 microloggers, which were located 
at 995, 1117, 1225, 1340, and 1485 m on the north face 
(Cailleret and Davi 2011), as follows (Table S2): 

r(z) = (PiTbasis + (p2^ + for temperatures (10) 

RH(z)= (Xi+X2-z) RHbasis+Zs for relative humidity 

(11) 

P(z) = (i//) + i/'2z)^'basis for precipitation (12) 

This 5-year climate sequence (from 2002 to 2006) was 
repeated in loops (six loops) for 30 years at each generation. 

Testing of hypotheses 

This study investigated how phenotypic plasticity and 
genetic adaptation contributed to TBB variation across ele- 
vations by comparing the following scenarios: 



Scenario A {'Neutral') was the baseline scenario without 
adaptive evolution, because -FcritBB was variable but not 
heritable {h^ = 0). Selection occurred within a genera- 
tion but was not expected to result in -FcritBB changes 
between generations. 

Scenario B ('Adaptive evolution), in which 

FcritBB was 

variable and heritable {h^ = 1) and selection occurred, 
potentially resulting in genetic evolution across genera- 
tions. 

Several variants of scenario B were also considered as fol- 
lows: 

Scenario C {'Evolution without mortality), in which 
FcritBB was variable and heritable {h^ = 1) and selection 
occurred only through differential reproduction without 
mortality. 

Scenario D {'Evolution without differential reproduction), 
in which FcritBB was variable and heritable {h^ = 1) and 
selection occurred only through differential mortality 
without differential reproduction among individuals. 
Scenario E {'Evolution, mortality driven by low level of 
cumulated carbon reserve'), in which mortality only 
occurred when the carbon reserve at the end of the year 
(CumCR) fell below a critical value (Type I mortality. 
Table 1). 

Scenario F {'Evolution, mortality driven by low level of car- 
bon reserve at budbursf), in which mortality only occurred 
when the carbon reserve before budburst (bbCR) feU 
below a critical value (Type II mortality, Table 1 ) . 
Scenario G {'Evolution, reduced heritability'), in which 
the heritability of FcritBB was set to h^ = 0.6. 
Scenario H {'Evolution with frost effect on LAP), in which 
every late frost reduced the LAI (Ha) by 10% or (Hb) by 
20% per degree below the critical minimal temperature. 

For each scenario, we ran 21 repetitions with different 
random initial conditions. Among the repetitions, only the 
spatial locations and the 500 initial founder genotypes chan- 
ged, whereas the allelic effects at each QTL were the same 
(Table 1). The average genetic value for FcritBB in the initial 
founder population was always = 190. The same alle- 

lic effects were used for all replicate runs of the simulation. 

Simulation result analyses 

The continuous population ranging between 700 and 
1700 m in elevation was split into five discrete adjacent 
populations, namely Altl (700-900 m), Alt2 (900- 
1100 m), Alt3 (1100-1300 m), Alt4 (1300-1500 m), and 
Alt 5 (1500-1700 m). Output variable distributions were 
obtained for each tree from the 30-year sequence (from 
ages 40 to 70). 

We computed the change in FcritBB and TBB between 
generations GO and G5 (Cb) population per population. As 
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^'critBB was Constant across the lifetime for a given tree, the 
change in -FcritsB was estimated as follows: 



(13) 



where n^ep is the number of repetitions (here n^^p = 5), and 
\^YnGx is the average -FcdtBB value at year n of generation X 
within the population under consideration (Y40 corre- 
sponded to the first year of the adult life stage). 

By contrast, TBB varied among the different climatic 
years (2002-2006). The change in TBB was estimated as 
follows: 



-Ef— E( 



l-'-y40+y,G5 l^y40+y,Go) 

(14) 



where «climYear = 5, and y = 0 for year 2002 up to y = 4 for 
year 2006. 

To measure the strength of within-generation selection, 
we also computed the change within each generation Gx as 

follows for FcritBB- 

Cwf„,b„ {Gx)=—J2{ 11f70Gx - liy40Gx ) (15) 

Note that both Cw and Cb are classically used in quanti- 
tative genetics; Cw is also referred to as the selection differ- 
ential (the difference of the mean trait value in a 
population before and after an episode of selection), and 
Cb measures the response to selection (the difference 
between the population distribution before selection and 
the distribution of the trait in the next generation). Cw and 
Cb values were compared between populations and scenar- 
ios using simple linear models with interaction between 
populations and scenarios (Appendix S4). All analyses were 
performed with R (RDevelopmentCoreTeam 2010). 

Results 

Because population Altl collapsed in almost all simulations 
(because of strong mortality below 800 m), it was excluded 
from the population-level results. 



Plastic response to the climatic gradient (scenario A) 

Only surviving trees were used for this analysis. Across all 
climatic years and all adult life stages, the length of growing 
season decreased on average from 210 to 160 days (Fig. 3A), 
and the WSI decreased by 33% between 1000 and 1700 m 
(Fig. 3B). As a consequence of these two limiting factors, 
the highest photosynthesis level (GPP = 1216 gc m^^) 
occurred at 1078 m (Fig. 3C) and was almost as high from 
-1050 to 1300 m. However, as respiratory costs also 



strongly decreased above 1200 m (Fig. 3D), the highest ring 
increment and seed production values were found at 
1258 m (Fig. 3E) and 1204 m (Fig. 3F), respectively. At ele- 
vations >1400 m, the ring width decrease was steeper than 
the seed production decrease. The minimal value of carbon 
reserves at the end of the year occurred between 1160 and 
1420 m, and the greatest difference between carbon reserves 
and carbon demand during budburst was found below 
1000 m (Fig. SI). Mortality was higher, on average, at low/ 
intermediate elevations (Table S3). 

The physiological response to the elevational gradient 
varied significantly among climatic years. For instance, ring 
widths were large regardless of elevation in 2002, but in 
2003, they increased continuously with the elevation, and 
in 2004, 2005, and 2006, they reached a maximum at 1000, 
1100, and 1200 m, respectively (Fig. S2A). This variability 
among years was also observed for seed production and the 
level of carbon reserves at the end of the year (Fig. S2B). 

Adaptive response of f critBB to the climatic gradient 

To investigate the effect of genetic adaptation along the 
gradient, we compared scenario B (adaptive evolution with 



1) with scenario A (neutral with hi 



0). First, 



we observed significant differences among scenarios and 
populations for the patterns of changes in -FcritBB between 
generations GO and G5 (denoted Cbp^^^gj), which measures 
the response to selection (Table 2, Fig. 4A,B, Appendix 
S4). In scenario A, Cbp._.„jj was low (<1°C) in all popula- 
tions. In scenario B, absolute Cbjj^^uj, values were signifi- 
cantly higher than in scenario A in populations Alt2 to Alt 
4 (Appendix S4). Populations Alt3 and Alt4 evolved toward 
lower -FcritBB values (on average, Cbp^^^jj = — 4.72°C in Alt3 
and -2.85°C in Alt4, Fig. 4), while population Alt2 evolved 
toward a higher _F„itBB (on average, Cbf^.^.^,„ = -l-1.18°C). 

Secondly, we investigated patterns of change in -FcritBB 
within each generation (Cw), as a measure of the strength 
of selection (Table 2, Fig. 5). Patterns of Cwp^^.^^j were sim- 
ilar among scenarios A and B at generation GO, with nega- 
tive Cwp^^„3j values within populations Alt3 to Alt5 
(Cwf^j-,ug = — 3.61°C, in Alt3, scenario B), and a positive 
Cwf^iujB value in population Alt2 (Cwp^^^jj, = +1.57°C, sce- 
nario B). Differences in absolute Cv/f^^.^^^ values indicate 
that selection was twofold more intense in population Alt3 
than Alt2 or Alt4, and it was weak in population AltS. 
These variations in selection direction and strength were 
consistent with the variations in FcntBB observed among 
generations (Cb) in scenario B. Finally, patterns of Cwf_.^„jb 
at generation G5 differed among scenarios A and B. While 
Cwf„„jB values remained similar across generation in sce- 
nario A, Cwp,^,|3j values decreased across generations in sce- 
nario B as a consequence of selection and recombination 
(Fig. 5). The phenotypic variances for FditBB within each 
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Figure 3 Elevational plastic variation in (A) lengtli of growing season (LGS), (B) water stress index (WSI), (C) gross primary production (GPP), (D) plant 
respiration (PR), (E) ring width (RW), and (F) seed production (SP). Each point corresponds to the average value across climatic years of the variable of 
interest for surviving trees (scenario B). 

Table 2. Simulated patterns of evolution for the temperature sum required for budburst (fcritBs)- For ^^ch population, the changes in FcntBB within 
generation GO (Cw) measures the intensity of selection, while the change between generations GO and GB (Cb) measures the response to selection. 
The phenotypic variance for Fcthbb (^'p) was computed at the last year of generation G5; in scenarios with = 1 (B to F; Ha and Hb), l/p is also the 
additive variance I/a. In scenario A, I/a = 0; in scenario G, Vp = I/a + Ve. Population Altl is not shown because of low population size. 



Population 

< 


Alt2 






Alt3 






Alt4 






Alt5 






Cb 


Cw 


l/p 


Cb 


Cw 


Vp 


Cb 


Cw 


l/p 


Cb 


Cw 


l/p 


A - Neutral 


-0.42 


1.32 


21.68 


0.48 


-4.28 


18.11 


-0.08 


-1.06 


21.13 


0.00 


-0.15 


21.37 


B - Adaptive evolution 


1.18 


1.57 


21.21 


-4.72 


-3.61 


13.60 


-2.85 


-1.11 


20.18 


-1.06 


-0.21 


19.64 


C - Evolution without mortality 


0.00 


0.00 


21.15 


-0.24 


0.00 


19.43 


-0.52 


0.00 


20.51 


-0.31 


0.00 


21.10 


D - Evolution without differential 


1.33 


1.36 


21.68 


-4.76 


-3.47 


14.55 


-2.69 


-0.81 


19.40 


-0.88 


-0.15 


20.41 


reproduction 


























E - Evolution, Type 1 mortality 


-0.23 


0.00 


19.64 


-4.93 


-3.72 


14.81 


-2.63 


-1.06 


19.03 


-0.94 


-0.18 


20.39 


F- Evolution, Type II mortality 


0.55 


0.98 


20.34 


-0.03 


-0.01 


19.81 


0.03 


0.00 


20.50 


0.06 


0.00 


19.98 


G - Evolution, reduced heritability 


0.60 


1.88 


33.72 


-3.79 


-6.09 


24.84 


-2.07 


-1.60 


32.31 


-0.90 


-0.35 


31.49 


Ha - Evolution, moderate effect of frost 


-1.21 


-0.40 


18.37 


-0.41 


-0.02 


20.30 


-3.71 


-2.07 


17.43 


-0.93 


-0.30 


18.37 


Hb - Evolution, strong effect of frost 


-1.15 


-0.39 


18.61 


-0.21 


0.03 


19.42 


-4.01 


-3.15 


18.21 


-1.17 


-0.33 


19.57 



population were also slightly lower in scenario B than in 
scenario A at generation G5, in particular in Alt3 (Table 2). 

Deciphering the mechanisms driving microevolution of 

fcritBB 

We investigated the effects of differential reproduction and 
differential mortality on evolutionary dynamics by compar- 
ing scenarios A, B, C, and D. Scenario without mortality 



(scenario C) led to weak changes in PcritsB from GO to G5 
(Fig. 4C), but overall the pattern of Cbp^.^,3j did not signifi- 
cantly differ from scenario A (Appendix S4), indicating 
that the absence of mortality prevented genetic adaptation. 
In contrast, scenario without differential reproduction (sce- 
nario D) resulted in a change in PcritBB from GO to G5 that 
was as important as it was in scenario B (Fig. 4D), indicat- 
ing that differential reproduction between trees played a 
minor role in simulated patterns of adaptation. 
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Figure 4 CInange in tine temperature sum required for budburst (Cb-FcniBB) between generations GO and G5 within populations Alt2 to AltS for dif- 
ferent scenarios (letter above each graph). The boxplot represents variation across the 21 repetitions. The dashed line correspond to Cb = 0 (no 
change). Population Altl was removed because of the low number of surviving individuals at G5. 



We also investigated the variations in Cbp^^,jj,j in scenario 
E, where mortality was driven by low levels of accumulated 
carbon reserves (Type I mortality), and scenario F, where 
mortality was driven by low carbon reserve levels at bud- 
burst (Type II mortality). In scenario E, populations Alt3 
and Alt4 still evolved significantly toward a lower f critsB as 
compared to scenario A, but the population Alt2 did not 
evolve toward higher faitsB (Fig- 4E, Appendix S4). This 
trend was caused by the absence of Type I mortality at low 
elevations (Table S3). In scenario F, population Alt2 still 
evolved toward a higher fcritBSi but populations Alt3 and 
Alt4 did not evolve toward a lower PcritBB (Fig. 4F). This 
result was due to the absence of Type II mortality in these 
populations (Table S3). The patterns of Cwp^^.^^^ values 
within populations in scenarios F and G were consistent 
with these Cb patterns (Fig. S3). 



Mortality not only drove evolution in fcritBB but also 
affected the elevational range of the population. In scenar- 
ios A and D, the range of the whole population (as mea- 
sured by the average elevation) shifted by on average -1-202 
and -1-168 m, respectively, between generations GO and G5 
(with minimal elevation -800 m), which could be due to 
higher mortality/lower reproduction at low elevation. By 
contrast, in scenarios C and E, the elevational shift of the 
whole population was reduced (-1-37 and +16 m respec- 
tively), and the minimal elevation at G5 was -720 m. 

Effect of heritability on genetic adaptation 

The heritability-level effects of fcritBB on the microevolu- 
tionary patterns of f critBB were analyzed by comparing sce- 
nario B {h^ = 1) to scenario G {h^ = 0.6). The divergence 
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Figure 5 Change in the temperature sum required for budburst (FcntBs) within generation GO (top panels) and GB (bottom panels), for scenarios A 
(neutral, left panels) and B (adaptive evolution, right panels). 



among populations that evolved positive (Alt2) versus neg- 
ative (Alt3 and Alt4) values for PcritsB was reduced when 
heritability was lower (Fig. 4). However, the response to 
selection remained important in populations Alt3 and Alt4 
(Cbf^,„3„ = -3.79°C and -2.07°C, respectively). 

Impact of frost on evolutionary dynamics 

We investigated the effects of frost using scenarios includ- 
ing a moderate (scenario Ha) or strong frost effect (sce- 
nario Hb) on the LAI. The patterns of microevolution 
markedly changed when compared to scenario B (Fig. 4). 
Population Alt2 no longer evolved toward higher values for 
fcritBB> and population 3 no longer evolved toward lower 
values of -FcritBB- Only population Alt4 evolved toward a 
lower F„i,BB value (Chp^.^^^ = -4.96°C and -3.55°C for 
scenarios Ha and Hb, respectively). The Cw patterns for 
J^critBB also Consistently indicated that selection was strong- 
est in population Alt4 (Fig. S3). Finally, frost effect pro- 
moted a reduced elevational shift of the whole population 
for scenarios Ha and Hb (+75.1 and +27.1 m, respectively) 



in comparison with scenario B (+167 m), which resulted 
from less severe mortality at lower elevations (Table S3). 

Discussion 

Effects of plasticity and microevolution on TBB variation 

Trees are long-lived species, and they experience a high var- 
iability in environmental conditions during their lives 
because of differences between the juvenile and adult stages 
and variations among climatic years. Accordingly, trees are 
expected to display a high plasticity for a wide range of 
functional traits, including TBB. In our baseline neutral 
scenario (A), the TBB varied by 12.2 days, on average, 
between the two extreme climatic years (and therefore 
7.6 days degree^') and by 35.2 days between the two 
extreme elevations (and therefore 5.4 days degree^'). This 
plastic variation was within the TBB range reported for 
_F. sylvatica by Vitasse et al. (2011), who found a range of 
TBB from 4.9 to 5.8 days degree^' across an elevational 
gradient from 131 and 1533 m in the Pyrenees. However, a 
lower plasticity was found in previous studies (between 2 
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and 2.5 days degree Kramer et al. 2008; Menzel et al. 
2001). 

In comparison with the plasticity, the TBB variation 
resulting from microevolution was small. Between the two 
most genetically differentiated populations (e.g., between 
populations Alt2 and Alt3 at generation G5) in the baseline 
scenario of adaptive evolution (B), the difference of 5.9°C 
in the temperature sum required for budburst (FcritsBi the 
genetic component of TBB) corresponded to an average of 
2 days for TBB (Fig. S4). These small variations are consis- 
tent with previous experimental studies; Vitasse et al. 
(2009b) found that a phenological model with constant 
parameters is able to reproduce TBB for different popula- 
tions, suggesting that plasticity hides local adaptation. 
Using a common garden experiment, Vitasse et al. (2009b) 
also reported that the genetic difference in TBB between 
populations originating from different elevations was 
almost 4 days, which is above our estimation of 2 days due 
to a wider elevational range. However, it is also likely that 
more differentiation would be found when simulating 
responses over more than five generations. 

Nonmonotonic elevation effects on the TBB genetic 
optimum 

Although they were weaker than phenotypic differentia- 
tion, significant patterns of genetic differentiation for TBB 
across elevations were nonetheless obtained in only five 
generations. In scenario B, one population (Alt2) evolved 
toward delayed budburst (-FcritBB increase of = -H.18°C), 
whereas populations Alt3 and Alt4 evolved toward earlier 
budburst (PcritBB decrease of -2.85°C to -4.72°C). The 
resulting pattern of -FcritBB variation across elevations was 
thus nonmonotonic. This result occurred because the mor- 
tality in population Alt2 was triggered by a lack of sufficient 
reserves before budburst to produce new leaves (Type II 
mortality, scenario E), and in populations Alt3 to Alt5, 
mortality was triggered by low reserves during the winter 
(Type I mortality, scenario F). The variability across eleva- 
tions in mortality-triggering factors resulted from the non- 
linear variability of the underlying ecophysiological 
processes, which makes the trees located at different eleva- 
tions pass different thresholds to mortality during different 
years. Type I mortality occurred mainly during 2004 for 
trees from populations Alt3 to AltS, and Type II mortality 
occurred mainly during 2003 for trees located below 
1100 m (i.e., populations Alt2 and Altl). In population 
Alt2, fewer trees displayed low reserves during winter in 
comparison with trees at higher elevations, because of 
longer growth period. However, these trees also displayed a 
higher LAI and biomass, which increased the respiratory 
and construction costs before budburst and made them 
more vulnerable to Type II mortality. 



When frost damage was considered (scenarios Ha and 
Hb), other nonmonotonic effects were observed. As only 
trees that initiated budburst could suffer from late frost, 
frost damage did not occur at the upper elevations (e.g., 
population Alt4) because delayed budburst protected them 
from late frosts. The first mechanistic consequence is that 
frost did not linearly affect the trees across the elevational 
gradient. Second, reduced LAI from frost damage can 
either increase or decrease the risk of mortality depending 
on its absolute effect for the reserve level. The risk of Type 
II mortality is expected to increase with frost damage 
because a reduced LAI decreases the carbon assimilation 
rate. However, a moderate frost in the model can also 
reduce Type I mortality because reduced LAI can decrease 
the carbon reserve required during budburst and the water 
loss during the following summer. This phenomenon 
explains why, when frost damage was considered, popula- 
tion Alt2 became less sensitive to carbon demand before 
budburst and did not evolve toward a later budburst as in 
the baseline scenario B. 

This study sheds light on the mechanisms that underlie 
genetic and phenotypic patterns of TBB variation. Consid- 
ering the number of underlying mechanisms involved in 
TBB and their patterns of environmental variation, this 
study suggests that nonmonotonic genetic patterns of TBB 
variation should be the rule rather than the exception. 
Other factors not considered here (for instance, assortative 
mating induced by variations in reproductive phenology 
across elevation) should be further studied (Soularue and 
Kremer 2012). 

Understanding the mechanisms underlying climate 
adaptation 

This study was based on a new mechanistic model coupling 
physiology, population dynamics, and quantitative genetics 
to simulate the short-term evolution of functional traits. 
Because PDG explicitly accounts for climate effects on the 
water and carbon exchanges of individual trees as a selective 
pressure, it provides a useful complement to existing evolu- 
tionary models of tree population life history traits (Le 
Corre and Kremer 2003; Kuparinen et al. 2010). More pre- 
cisely, in PDG as in these other models, individual fitness is 
the parameter driving the process of adaptation. However, 
the primary originality of PDG is that individual fitness is 
an output, which is calculated as the lifetime reproductive 
success resulting from a combination of functional traits 
and the environmental context. This was also used (Kramer 
et al. 2008) to investigate the temporal patterns of micro- 
evolution in a single population. We extended this 
approach here, and we showed how such a mechanistic 
model can be used to investigate the type and strength of 
selection mediated by the climate through the estimation of 
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the selection differential (using Cw, the difference in the 
population before and after an episode of selection). 

In PDG, selection occurred both through differential 
mortality and reproduction of individual trees within each 
generation. Mortality was found to be the main driver of 
evolutionary dynamics, with different types of mortality 
promoting different patterns of adaptation across eleva- 
tions. Predicting tree mortality is a key and complex issue 
in tree physiology and ecology because many mechanisms 
are involved and interact (carbon starvation, cavitation, 
and pathogens, (McDowell et al. 2008). We chose to model 
tree mortality according to the carbon starvation hypothe- 
sis alone, in which a tree dies when carbon reserves are too 
low to allow the setup of new leaves in the spring (second 
threshold. Type II mortality) or to ensure tree functioning 
during the winter (first threshold. Type I mortality). We 
excluded two other mechanisms because (i) no major 
pathogens were observed in our site for F. sylvatica and (ii) 
in Ventoux the minimal hydraulic potential (—2 Mpa) is 
above the critical pressure causing a 50% loss of conduc- 
tance (-2.4 Mpa, Herbette et al. 2010). 

In contrast, differential reproduction among individuals 
was found to have a minor role as a driver of evolutionary 
dynamics. To our knowledge, few models relate reproduc- 
tion to tree carbon cycles (Genard et al. 2008). Our main 
hypothesis here was that seed production increases with car- 
bohydrate reserves. This idea is consistent with the higher 
seed production observed for dominant trees (H. Davi, per- 
sonal observation) and for years after good climatic years. It 
is also consistent with the resource supply hypothesis, in 
which fi-uit production, especially during mast years, occurs 
when carbohydrate reserves are sufficient (Yamauchi 1996). 

In addition to selection, the outcome of adaptation is also 
determined by the level of genetic variation available for 
selection, which depends on the heritability of the trait under 
selection and on genetic architecture (in particular, the num- 
ber of QTLs determining the trait's genetic variation). We 
showed that levels of heritability for the TBB such as those 
measured for F. sylvatica (h^ = 0.6, Kramer et al. 2008) lead 
to significant patterns of genetic differentiation for the TBB 
across elevations. It was out of the scope of this study to 
investigate genetic architecture effects in detail. Therefore, 
we chose a simple additive quantitative genetic model with 
ten independent QTLs for TBB, as in previous studies (Le 
Corre and Kremer 2003; Kuparinen et al. 2010). However, 
further investigations into the effects of the quantitative 
genetic model are needed (Le Corre and Kremer 2012). 

Main shortcomings of the mechanistic Physio-Demo- 
Genetic model 

Admittedly, fuU PDG evaluation requires the parameteriza- 
tion and evaluation of complex mechanisms of tree physi- 



ology, demography, and selection that was beyond the 
scope of this study. Among the most important shortcom- 
ings of this study, we repeatedly reused a short (5-year) cli- 
matic period that is clearly unrealistic and potentially 
biases the estimates of survival between generations. More- 
over, this period included specific climatic years, which can 
have major influence on results. For instance, 2004 was an 
exceptional drought year that has led to low growth rates of 
beeches (Cailleret and Davi 2011). Second, PDG does not 
include competition between adult trees, which is a process 
that is potentially more important for tree growth and sur- 
vival than climate. However, not accounting for competi- 
tion in studying microevolution driven by TBB is 
reasonable because TBB is related primarily to temperature 
and only secondarily to competition. Indeed, we previously 
showed that dominant trees exhibited an earlier TBB, but 
this effect is small in comparison with the elevation and 
year effects (Davi et al. 2011). Third, we used nonoverlap- 
ping generations, which make individual level comparisons 
with forest inventory data and tree ring increment mea- 
surements difficult. Fourth, dormancy was not taken into 
account in modeling F. sylvatica budburst. Fifth, a real sen- 
sitivity analysis on the entire model (as performed for the 
ecophysiological component of PDG in Dufrene et al. 
2005) will be needed to strengthen some of our conclusions 
and to draw a more accurate picture of what process con- 
trol genetic adaptation of budburst. 

Nevertheless, valid conclusions could be drawn using the 
current version of PDG. This success is possible mainly 
because the physiological module of PDG (CASTANEA), 
which models the climate effect on tree functioning, has 
already been thoroughly validated for European beech in 
several previous studies (e.g., Dufrene et al. 2005). PDG 
simulated the maximum tree ring width for elevations 
between 1100 and 1420 m, which also corresponds to the 
maximum ring width observed in the field (Cailleret and 
Davi 2011). 

Conclusions 

We described here a new modeling tool (PDG) to assess 
the potential mechanisms of local adaptation for trees 
under changing environmental conditions. The primary 
originalities of the PDG model are that it combines physi- 
ology, demography, and genetics and that fitness is a 
dynamic output of the model. Such complex models are 
useful tools for predicting the evolution of nonequilibrium 
forest populations under CC, under which many tipping 
points and nonlinear effects may be involved. PDG model 
requires a large amount of data to be parameterized and 
tested, and results must be cautiously interpreted. Demo- 
graphic processes such as mortality and reproduction 
should be further studied, and other processes such as 
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competition and regeneration have to be included in this 
general framework. Nevertheless, the following two impor- 
tant conclusions emerge from our present study: (i) 
Genetic evolution of tree populations can occur in a few 
generations (<5), and (ii) patterns of genetic differentiation 
across space (and across elevations here) can be nonmono- 
tonic. 
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